{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Construction and meshing of fracture networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial, we will show:\n",
    "\n",
    "1. How to define fractures and a fracture network in a 3d domain.\n",
    "2. How to construct a family of meshes that represent the 3d domain, the fractures and their intersections.\n",
    "3. Assembly of the grids into a `GridBucket` container that stores all grids, and the geometric relation between them.\n",
    "\n",
    "Together, these are the first steps towards creating a simulation model for a mixed-dimensional problem in fractured domains.\n",
    "\n",
    "\n",
    "## Summary\n",
    "For most simulation purposes, the final grid bucket is all that is needed. Therefore, we start by showing a shortcut for obtaining a `GridBucket` given a set of fractures, a domain and two mesh size parameters. All these will be described in more detail below.\n",
    "\n",
    "### Meshing of 2d fractures\n",
    "Fracture networks in 2d are formed by a set of points, together with their connections. We define this by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import porepy as pp\n",
    "\n",
    "# Point coordinates, as a 2xn array\n",
    "p = np.array([[0, 2, 1, 1], [0, 0, 0, 1]])\n",
    "# Point connections as a 2 x num_frac arary\n",
    "e = np.array([[0, 2], [1, 3]])\n",
    "\n",
    "# The domain contains two fractures: The first from (0, 0) - (0, 2), the second (1, 0) to (1, 1)\n",
    "# Set domain boundaries\n",
    "domain = {'xmin': -2, 'xmax': 3, 'ymin': -2, 'ymax': 3}\n",
    "\n",
    "# Define a fracture network in 2d\n",
    "network_2d = pp.FractureNetwork2d(p, e, domain)\n",
    "\n",
    "# Set preferred mesh size close to the fracture, and at the boundary (essentially this is a far-field value)\n",
    "mesh_args = {'mesh_size_frac': 0.2, 'mesh_size_bound': 0.3}\n",
    "\n",
    "# Generate a mixed-dimensional mesh\n",
    "gb = network_2d.mesh(mesh_args)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3d meshing\n",
    "Fractures in 3d are polygons instead of lines. This makes computations of the geometry quite a bit more difficult, and requires further data structure to store the fractures:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The fractures are specified by their vertices, stored in a numpy array\n",
    "f_1 = pp.Fracture(np.array([[0, 1, 2, 0], [0, 0, 1, 1], [0, 0, 1, 1]]))\n",
    "f_2 = pp.Fracture(np.array([[0.5, 0.5, 0.5, 0.5], [-1, 2, 2, -1], [-1, -1, 2, 2]]))\n",
    "\n",
    "# Also define the domain\n",
    "domain = {'xmin': -2, 'xmax': 3, 'ymin': -2, 'ymax': 3, 'zmin': -3, 'zmax': 3}\n",
    "\n",
    "# Define a 3d FractureNetwork, similar to the 2d one\n",
    "network = pp.FractureNetwork3d([f_1, f_2], domain=domain)\n",
    "mesh_args = {'mesh_size_frac': 0.3, 'mesh_size_min': 0.2}\n",
    "\n",
    "# Generate the mixed-dimensional mesh\n",
    "gb = network.mesh(mesh_args, ensure_matching_face_cell=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the option ensure_matching_face_cell=False - this is needed to circumvent some issues with the way gmsh produces grids. As far as we are aware, it has no consequences for the simulation quality, but we are working on fixing it.\n",
    " \n",
    "## Import filters \n",
    "FractureNetworks (2d and 3d) can also be defined directly from files storing their data, see pp.fracture_importer for details. \n",
    " \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Technical details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We focus on the the technicalities related to 3d meshing; in 2d, meshing is relatively simple (still difficult, but 3d is much worse). \n",
    "\n",
    "Functionality for fractures and their intersection are provided in the subpackage `porepy.fracs`. Fractures are defined either as Elliptic fractures, or as convex, planar polygons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify a fracture by its vertexes, as a 3xn array\n",
    "pts_1 = np.array([[0, 1, 2, 0], [0, 0, 1, 1], [0, 0, 1, 1]])\n",
    "f_1 = pp.Fracture(pts_1)\n",
    "\n",
    "# .. and another fracture, intersecting the first\n",
    "pts_2 = np.array([[0.5, 0.5, 0.5, 0.5], [-1, 2, 2, -1], [-1, -1, 2, 2]])\n",
    "f_2 = pp.Fracture(pts_2)\n",
    "               "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also specify the fracture as an ellipsis, approximated as a polygon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify the fracture center\n",
    "center = np.array([0.1, 0.3, 0.2])\n",
    "# The minor and major axis\n",
    "major_axis = 1.5\n",
    "minor_axis = 0.5\n",
    "\n",
    "# Rotate the major axis around the center.\n",
    "# Note that the angle is measured in radians\n",
    "major_axis_angle = np.pi/6\n",
    "\n",
    "# So far, the fracture is located in the xy-plane. To define the incline, specify the strike angle, and the dip angle.\n",
    "# Note that the dip rotation is carried out after the major_axis rotation (recall rotations are non-commutative).\n",
    "strike_angle = -np.pi/3\n",
    "dip_angle = -np.pi/3\n",
    "\n",
    "# Finally, the number of points used to approximate the ellipsis. \n",
    "# This is the only optional parameter; if not specified, 16 points will be used.\n",
    "num_pt = 12\n",
    "f_3 = pp.EllipticFracture(center, major_axis, minor_axis, major_axis_angle, strike_angle, dip_angle, num_points=num_pt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fractures can be joined into a `FractureNetwork`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "network = pp.FractureNetwork3d([f_1, f_2, f_3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The `FractureNetwork` class is the base for analysis and manipulation of fracture networks. The functionality is expanding on demand. For the moment, the most interesting feature is the export of the fracture network to ParaView (requires the vtk extension of python installed, see installation instruction):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/eke001/anaconda3/envs/porepy/lib/python3.6/site-packages/vtk/util/numpy_support.py:135: FutureWarning: Conversion of the second argument of issubdtype from `complex` to `np.complexfloating` is deprecated. In future, it will be treated as `np.complex128 == np.dtype(complex).type`.\n",
      "  assert not numpy.issubdtype(z.dtype, complex), \\\n"
     ]
    }
   ],
   "source": [
    "network.to_vtk('fracture_network.vtu')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting file can be opened in ParaView. A little bit of work in ParaView gives the following picture"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"fracture_network.png\" width=300>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have not yet set a boundary for the `FractureNetwork`, and effectively for the domain. The boundary is defined as a box, and is imposed in the following way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The domain is a dictionary with fields xmin, xmax, etc.\n",
    "domain = {'xmin': -2, 'xmax': 3, 'ymin': -2, 'ymax': 3, 'zmin': -3, 'zmax': 3}\n",
    "network.impose_external_boundary(domain)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above, we defined the bounding box to not intersect with the fractures. If the domain would have been smaller, fractures that intersect a face of the box would by default (can be overruled) have been truncated so that they are confined within the bounding box."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Meshing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our aim is to create a computational mesh that conforms to the fractures, as well as to their intersections (1d lines, 0d points). For the actual grid construction we rely on Gmsh. However, these packages all require that the geometric constraints, that is the fractures, are described as *non-intersecting* polygons [if you know of packages that do not require this, please let us know]. It only takes some thinking to understand why the meshing software would not like to do this themselves; this is a highly challenging task.\n",
    "\n",
    "PorePy provides functionality for finding intersections between fractures, and splitting them into non-intersecting polygons. Intersections are found by "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "network.find_intersections()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get information on the number of intersections, type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'In total 9 fractures intersect in 15 intersections'"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "network.intersection_info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fracture and intersection count also includes the bounding planes of the domain, hence the number 9 instead of 3.\n",
    "\n",
    "When we have found all intersections, the fracture planes should be split into polygons that do not intersect, but that may share edges along intersection lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "network.split_intersections()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geometric tolerances and stability of meshing algorithm\n",
    "A critical concept in meshing of fractured domains is the concept of geometric tolerance: Since we are operating in finite precision aritmethics, two points may or may not be consider equal (or similarly, two lines / planes may or may not intersect), depending on how accurately we consider their representation. At least three concepts come into play here\n",
    "\n",
    "1. The accuracy of the numerical representation of the objects (accumulated effect of finite precision rounding errors).\n",
    "2. The accuracy in geological interpretation of fracture data: If the fracture network originates from an interpretation of satellite images, differences measured in centimeters should be treated with some caution\n",
    "3. The resolution of the computational grid: If points with a certain distance are considered non-equal, this may also require that we resolve their difference in the mesh. In addition, the mesh generator will use its own concept of geometric tolerance for internal calculations.\n",
    "\n",
    "In PorePy, these issues are attempted resolved as follows: The `FractureNetwork` has an attribute `tol` that represent the geometric tolerance used in the calculation of intersections and subsequent splitting of the fractures. If meshing is done with gmsh, the tolerances used in PorePy and gmsh are related. The approach works reasonably well, but for complex configurations of fracture intersections, stability issues can arise. We are working to iprove these matters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interaction with gmsh\n",
    "\n",
    "Next, create grids for the domain, as well as for fractures and fracture intersections. This involves creating a config file for the mesh generator that contains geometry description, including fracture planes and their intersections. The mesh is then created by calling gmsh (NOTE: The path to the gmsh executable should be specified in a PorePy config file, type 'porepy.utils.read_config?' for more information). The resuling mesh information is read back to python, and `Grid` objects representing the matrix, fractures and fracture intersections are created.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gmsh is quite flexible in terms of letting the user set / guide the preferred mesh size in different parts of the domain. PorePy tries to adjust to this adapting the specified mesh size to the geometry. From the user side, two parameters must be specified: mesh_size_frac gives the target mesh size in the absence of geometric constraints, while mesh_size_min gives the minimal mesh size to be specified to Gmsh. What actually happens with the mesh, that is, how Gmsh translates these preferred options into a grid, is another matter. It may take some practice to get this to work properly.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of cells: 25938\n",
      "Number of faces: 55023\n",
      "Number of nodes: 5819\n"
     ]
    }
   ],
   "source": [
    "mesh_size_frac = 0.3\n",
    "mesh_size_min = 0.2\n",
    "mesh_args = {'mesh_size_frac': mesh_size_frac, 'mesh_size_min': mesh_size_min}\n",
    "\n",
    "# With the mesh size parameters, we can simply ask the GridBucket to mesh itself:\n",
    "gb = network.mesh(mesh_args)\n",
    "\n",
    "# Report the number of cells, faces and nodes before and after.\n",
    "g = [g for g, _ in gb if g.dim == gb.dim_max()][0]\n",
    "print('Number of cells: ' + str(g.num_cells))\n",
    "print('Number of faces: ' + str(g.num_faces))\n",
    "print('Number of nodes: ' + str(g.num_nodes))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The GridBucket is in effect a mixed-dimensional mesh that can be used for discretization and visualization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualization of the mixed-dimensional mesh\n",
    "The set of meshes in the `GridBucket` can be dumped to ParaView by simply writing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/rbe051/anaconda3/envs/porepy/lib/python3.6/site-packages/vtk/util/numpy_support.py:134: FutureWarning: Conversion of the second argument of issubdtype from `complex` to `np.complexfloating` is deprecated. In future, it will be treated as `np.complex128 == np.dtype(complex).type`.\n",
      "  assert not numpy.issubdtype(z.dtype, complex), \\\n"
     ]
    }
   ],
   "source": [
    "e = pp.Exporter(gb, 'grid_bucket')\n",
    "e.write_vtk()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, some manipulations in ParaView show how the grids on fracture surfaces intersects with the matrix grid.\n",
    "\n",
    "<img src='mixed_dimensional_grid.png'  width=200>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Next steps\n",
    "Now that we have created the `GridBucket`, the next step is to solve mixed-dimensional flow and transport problems. This is covered by the tutorial Darcy_equation, among others."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
